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Abstract. We derive an estimate of the rate of formation of dark matter haloes per unit 
volume as a function of the halo mass and redshift of formation. Analytical estimates 
of the number density of dark matter haloes are useful in modeling several cosmolo- 
gical phenomena. We use the excursion set formalism for computing the formation rate 
of dark matter haloes. We use an approach that allows us to differentiate between major 
and minor mergers, as this is a pertinent issue for semi-analytic models of galaxy form- 
ation. We compute the formation rate for the Press-Schechter and the Sheth-Tormen 
mass function. We show that the formation rate computed in this manner is positive at 
all scales. We comment on the Sasaki formalism where negative halo formation rates 
are obtained. Our estimates compare very well with N-body simulations for a variety 
of models. We also discuss the halo survival probability and the formation redshift 
distributions using our method. 

Keywords : cosmology: large scale structure of the Universe - galaxies: formation - 
galaxies: evolution - galaxies: haloes 

1. Introduction 



Gravitational amplification of density pe rturbations is thought to be responsible for formation 
of large scale structures in the Universe ([Peebles Il980[ IShandarin & Zeldovichl l Il989t iPeacockl 
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1999 : Padmanabhanl 2002 : Bemardeau et al. I I2OO2I ). Much of the matter is the so called dark 



matte r that is believed to be weakly interacting and non-relativistic (Trimble 1987; K omatsu et al 



20091). Dark matter responds mainly to gravitational forces, and by virtue of a larger density than 



baryonic matter, assembly of matter into haloes and large scale structure is primarily driven by 
gravitational instability of initial perturbations in dark matter Galaxies are believed to for m when 
gas in highly over-dense haloes cools and collapses to form stars in significant numbers dHovlel 



, Lap; _ 

I953I: iRees & Ostrike3ll977t IsilSl 19771: lBinnevl[l977l) . Thus the hierarchical formation of dark 



matter haloes is the key driver that leads to formation and evolution of galaxies and clusters of 
galaxies. 



The halo mass function describes the comoving number density of dark matter haloes as 
a function of mass and redshift in a given cosmology. It is possible to develop the theory of 
mass functions in a manner that makes no reference to the details of the cosmological model 
or the power spectrum of fluctuations. That is, we expect the mass function to take a univer- 
sal form, when scaled appropriately. Simple theor etical arguments have been used to obtain 
this uni versal functional form of the mass function feess & Schechter"l974'; 'Bond et aL 1991 ; 
Sheth. M o & Tormen 2001). Bond et al. ( 1993), andTIsheth, Mo & Tormen (2001) used the ex- 



cursion set theory to derive the mass function. Much work has also been done to determine 
the extent to which this form is consistent with results from N-body simulations ("Jenkins et al 
2001; White 2002; Reed et al. 2003; W arren et al. 2006; Reed e t al. 2007; Lukic et al. 2o6" 



lllC 



Cohn & Whitell2008t iTinker et al.ll2008l) with the conclusion that the agreement is fairly good. 



It is remarkable that a purely local approach provides a fairly accurate description of the mani- 
festly non-linear and strongly coupled process of gravitational clustering. The success of the 
local description has been exploited in developing the semi-analytic theories of galaxy formation 



jWhite & Frenk''l99l';'Ka uffmann. White & Guiderdonill993l;IChiu & Ostrikej2000[|Madau. Ferrara & Rees 
|20pi;Samui et al. 2007). 



The Press-Schechter mass function (I Press & Schechterl 19741) that is common ly used in these 
semi-analytic models assumes spherical collapse of haloes dOunn & Gott 1972 ). The shape of 
this mass fun ction agrees with num e rical results qualitat ively, but there are deviations at a quant- 



iilts qu alitati' 
et al.l l200lV 



itative level (lEfstathiou et al.lll988t iJenkins et al.l 12001). Improvements to the Press-Schechter 
mass function have been made to overcome this limitation. In particular, t he Sheth-Tormen 



mass function, which is based on the more realistic ellipsoidal collapse model ( Sheth & TormenI 
19991; ISheth. Mo & Tormenll200lb fits numerical results better Many fitting functions with three 



or four fitting parameters have been propose d, these are based on results of simulations of the 



Lambda-Cold Dark Matter (ACDM) niodel (IJenkins et al.ll200U iReed et alJl2003l: IWarren et al 
20061: iFakhouri. Ma & Bovlan-KolchirJfoiol) . 



In the application of the theory of mass functions to the semi-analytic models for galaxy 
formation, we often need to know comoving number density of haloes of a certain age. Naturally, 
this quantity is related to the halo formation rates and the survival prob ability. While these detail s 
are known and well understood for the Press-Schechter mass function (iPress & Schechtedl 19741) . 
the situation is not as clear for other models of the mass function. Furthermore, analytic estimates 
for the halo formation rate and survival probability are important in spite of the availability of 
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accurate fitting functions for these quantities in the ACDM model. This is because analytic 
estimates can be used to study variation in these quantities with respect to, for instance, the 
underlying cosmology or the power spectrum of matter perturbations. Studying such variation 
with the help of simulations is often impractical. In this work, we focus on the computation of 
halo formation rates. 



Several approach es to calculating halo formation rates have be en sug gested (IB lain & Longaii 
199llSasak i 1994; Ki tavama & Sutoll 19961) . In particular. ISasakil ( 1 19941) suggested a very simple 



approximation for the formation rate as well as survival probability for haloes. The approximation 
was suggested for the Press-Schechter mass function, though it does not use any specific aspect 
of the form of mass function. The series of arguments is as follows: 



• Merger and accretion lead to an increase in mass of individual haloes. Formation of haloes 
of a given mass from lower mass haloes leads to an increase in the number density, whereas 
destruction refers to haloes moving to a higher mass range. The net change in number 
density of haloes in a given interval in mass is given by the difference between the formation 
and destruction rate. 

• Given the net rate of change, we can find the formation rate if we know the destruction rate. 

• A simple but viable expression for the destruction rate is obtained by assuming that the 
probability of destruction per unit mass (also known as the halo destruction efficiency) is 
independent of mass. 

• This approximate expression for the destruction rate is then used to derive the formation 
rate as well as the survival probability. 



The resulting formulae have been applied freely to various cosmologies and power spec- 
tra, including the CDM class of power sp ectra. The S asaki approa ch has been used in many 



semi-analytic mod els of galaxy formation (iChiu & Os triker 2000; C houdhurv & Ferraral 12005 



Samui et alJl2007n mainly due to its simplicity. Attempts have also been made to generalize 
the approxirnation to models of mass function other than the Press-Schechter mass function 
dSamui et al ]l2009l) . though it has been found that a simple extension of the approximation some- 



times leads to unphysical results. In particular, when applied to the Sheth-Tormen mass function, 
the Sasaki approach yields negative halo formation rates. 

In this paper, we investigate the application of the Sasaki approach to the Sheth-Tormen mass 
function. We test the Sasaki approach by explicitly computing the halo formation and destruction 
rates for the Press-Schechter mass function using the excursion set formalism. We then generalize 
this same method to compute the halo formation rates for the Sheth-Tormen mass function. We 
find that halo formation rates computed in this manner are always positive. Finally, we compare 
our analytical results with N-body simulations. 



A reason for choosing the approach presented in this paper, as compared to other compet- 
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ing approaches based on the excursion set formalism, is that we wish to be able to differen- 
tiate between major and minor mergers. This is an essential requirement in semi-analytical 
models of galaxy forma t ion and is not addressed by other approaches for halo formation rate 
jPercival & Miller" 1999'; 'Percival. Miller & PeacocM l2000t IPercivall I2OO1I lOiocoli etal.1 12007 
iMoreno. Giocoli & Sheth 2008, 200^^ 



Many pre vious studies of merg er rates using analytical or numerical techniques are present in 
the literature. iBenson et al.l (l2005b recognised that the Sasaki approach of calculating halo form- 
ation rate was fundamentally inconsistent. They showed that a mathematically consistent halo 
merger rate should yield current halo abundances when inserted in the Smoluchowski coagulation 
equation. They applied this technique to obtain merger rates for the Press-Schecht er mass func- 
tion. The original formulati on of halo merger r ates in the excursion set picture (jL acey & Cole! 
1994 was also improved by iNeistein & Dekell (l2008b and iNeistein. Maccio & Dekel C201Q) to 
include the effect of finite merger time inter val. They found that the resultant merger rates are 
about 20% more accurate than the estimate of iLacev & Coi3 (Il994 for minor mergers and about 
three times more accura te for minor me rgers. However, most of these studies have focussed on 
the overall merger rates ( Cohn & Whitel l2008) rather than halo formation rates. 



We discuss the Sasaki and the excursion set formalisms in §2. We describe our simulations 
in §3, discuss our results in §4 and finally summarise our conclusions in §5. 



2. Rate of halo formation 



The total change in number density of collapsed haloes at time t with mass between M and M+dM 
per unit time is denoted by N{M, t)dM and is due to haloes gaining mass through accretion or 
mergers. Lower-mass haloes gain mass so that their mass is now between M and M + dM, and 
some of the haloes with mass originally between M and M + dM gain mass so that their mass 
now becomes higher than this range. We call the former process halo formation and the latter 
as halo destruction, even though the underlying physical process is the same in both cases; the 
different labels of formation or destruction arise due to our perspective from a particular range of 
mass. We denote the rate of halo formation by A^form(A/, t)dM and the rate of halo destruction by 
NAe!,t(M, t)dM. We immediately have 

N{M, t) = A?form(M, f) - A?dest(M, f)- (1) 



Following Sasakil ( 1994 ). in general we can formulate each term in the above expression as fol- 



lows. The rate of halo destruction can be written as 



A^dest(M,f) = N(M,t)Q(M,M';t)dM' (2) 



M 

= (f>{M, t)N(M, t), (3) 

where, Q{M, M'; f) represents the probability of a halo of mass M merging with another halo to 
form a new halo of mass M' per unit time. The fraction of haloes that are destroyed per unit time 
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is denoted by 4>{M, t). This quantity is also referred to as the efficiency of halo destruction. The 
rate of halo formation can be written as 

M 

NformiM, t)^ j N(M' , t) Q{M' , M; t)dM' (4) 



where Q{M',M; t) represents the probability of a halo of mass M' evolving into another halo of 
mass M per unit time. We can now write, from equation ([T} and from our definitions in equations 
© and ©, 

A^form(M, r) = NiM, t) + 0(M, t)N{M, t). (5) 
This reduces the calculation of rate of halo formation to a computation of 4>{M, t). 



Sasakil (119941) proposed a simple ansatz to compute (f>{M, t): if we assume that the efficiency 



of halo destruction has no characteristic mass scale and we require that the destruction rate re- 
mains finite at all masses then it can be shown that (p does not depend on mass. 



2.1 Sasaki prescription: Press-Schechter mass function 



To understand the Press-Schechter formalism dPress & SchechteJ 19741: iBond et "Zl ll99lh . which 



gives the co-moving number density of collapsed haloes at a time f with mass between M and 
M + dM, let us consider a dark matter inhomogeneity centered around some point in the Universe. 
The smoothed density contrast within a smoothing scale of radius R around this point is defined 
as S(R) = [p(R) - p]/p, where p{R) is the density of dark matter within R and p is the mean 
background density of the Universe. If this density contrast 6{R) is g reater than the thres hold 
density contrast for collapse 6c obtained from spherical collapse model ( Gunn & Gotdll972 ). the 



matter enclosed within the volume collapses to form a bound object. In hierarchical models, 
density fluctuations are larger at small scales so with decreasing R, 6{R) will eventually reach 6c- 
The problem then is to compute the probability that the first up-crossing of the barrier at 6c occurs 
on a scale R. This problem can be addressed by excursion set approach. 

The excursion set approach consists of the following principles: consider a trajectory 6{R) 
as a function of the filtering radius R at any given point and then determine the largest R at 
which 6{R) up crosses the threshold 6c{t) corresponding to the formatio n time t. The solutio n 



of the problem can be enormously simplified for Brownian trajectories (IChandrasekharlll943l) . 
that is for sharp fc-space filtered density fields, as in this case contribution of each wave mode is 
independent of all others. In such a case we have to solve the Fokker-Planck equation for the 
probability density Tl(6,S)d6, where S = cr^{R) and cr{R) is the standard deviation of fluctuations 
in the initial density field, smoothed at a scale R, 

dn{6,S) _ 1 d^n{6,S) 

as " 2 86^ ^ ' 
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The solution (IPorciani et al.lll998 ; IZentneilEoOTi) can be obtained using the absorbing boundary 
condition \lii6c{t),S) = and the initial condition Il{5, S - 0) - 6d(S), where SoiS) is the Dirac 
delta function 



n(6, S ; 6c)d6 



1 



JIjtS 



exp|-— |-exp| — 



IS 



d6 



(7) 



Now define F{S, 6c{t)) - j d6Tl(6, S ; 5c(t)) as the survival probability of trajectories and ob- 
tain the differential probability for a first barrier crossing: 



f(S) 



dF(S,6At)) 6At) 



dS 



: exp 



2S 



(8) 



From this, one can obtain the co-moving number density of collapsed haloes at time t with mass 
between M and M + dM 



NpsiM,t)dM = ^f(S) 



dS 



dM 



dM 



2Pm- , 1 

"77 W 

71 M 



d In cr 




V 


dM 


exp 


~2_ 



dM. 



(9) 



here p^r is the average comoving density of non-relativistic matter and v = [Scit)/o-{M))]^ = 
[6c/(D{t)o-{M))]^, where 6c is the threshold density contrast for collapse, D{t) is the linear rate of 
growth for density perturbations and cr(M)(= S '^^) is the standard deviation of fluctuations in the 
initial density field, which is smoothed over a scale that encloses mass M. 

In the following discussion, we will denote the mass function by N(M, t) if the statement 
is independent of the specific form of the mass function. We will use a subscript PS when the 
statements apply only to the Press-Schechter form of the mass function. 



With Sasaki's ansatz, the destruction rate efficiency (f) can be written in this case as 
1 dD{t) 



Dif) dt 



(10) 



With this, we can write down the rate of halo formation for the Press-Schechter mass function 
from equation (|5]i as: 



1 dDit) 



1 dD{f) 
IXt) dt 



D(t) dt 
Nps{M,t)- 



^(M)DHty 



(11) 



Note that for haloes with large mass, that is in the limit 6c s> cr(M)D{t), Nform approaches A^ps. In 
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other words, the total change in the number of haloes is determined by formation of new haloes. 
For haloes with low mass, where cr is much larger than unity, although A^form remains positive, 
the total change is dominated by destruction and A^ps becomes negative. 

We can also define two related, useful quantities now. Firstly, the probability p(t\ , t2) that a 
halo which exists at f i continues to exist at f2 without merging is given by 



p(ti,t2) = exp 



'2 

■ J cf>(t')dt' 



D(h) 

D(t2) 



(where t2 > fi) 



(12) 



This is usually known as the survival probability of haloes, and is independent of halo mass in 
the Sasaki prescription. In this picture, the distribution of epochs tf of formation of haloes with 
mass M at time t is given by 



F(M; tf, t)dMdtf - NfoimiM, tf)p{tf, t)dMdtf. 



(13) 



2.2 Sasaki prescription: Sheth-Tormen mass function 



The Press-Schechter mass function does not provide a very good fit to halo mass function obtained 
in N-body simulations. In particular, it under-predicts the number density of large mass haloes, 
and over-predicts that of small mass haloes. Hence it is important to generalize the calculation 
of formation rates to other models for mass function that are k nown to fit simulations better The 
Sheth-Tormen form of mass function ( Sheth & TormenI 1999 ) is known to fit simulations much 
better than the Press-Schechter form£ (For a comparison of both of these forms of halo mass 
function with simulations, see Fig. 3 of iJenkins et al ] |200I[ .) The Sheth-Tormen mass function is 
given by 



A^st(M, t)dM ^AJ-— (avy^ 

71 M 



2pn 



d In cr 



dM 

X [1 + (avy''] exp 



av 



dM, 



(14) 



where the parameters a , p, and A have best fit values of a = 0.707, p = 0.3 and A = 0.322 
dSheth & Tormenlll999h . and the quantity v is as defined before. This form of mass function has 
the added advantage of being similar to the mass function derivedjising a varia ble barrier motiv- 
ated b y ellipsoidal collapse of overdense regions dSheth. Mo & TormenI boOll ISheth & TormenI 
20021) . Note that if we choose A = 0.5, p = and a - 1 then we recover the Press-Schechter mass 



Even this fomi of halo mass function has poor accuracy in some cases, namely, for conditional mass functions with 
large mass ratios and for mass function in overdense regions (Sheth & Tormen 2002). In applications involving these 
regimes it is perhaps advisable to use more accurate fitting functions to simulation data. However, the Sheth-Tormen form 
still has the property of being considerably better than the Press-Schechter form while having a physical interpretation. It 
is thus preferable in many semi-analytic models where the Press-Schechter form is used. 
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function derived using spherical collapse. Recently, it has been shown that the best fit values o f 
these parameters depend on the slope of the power spectrum (Bagla, Khandai & Kulkarnil 20091) . 

We can now apply the Sasaki p rescription to this form of mass function and calculate the rates 
of halo formation and destruction ( Ripamontill2007 1. We get for the destruction rate efficiency 



^« = 7^^[l-2/']- (15) 
D at 

Note that the destruction rate efficiency is independent of mass. The rate of halo formation is then 
given by 



■ ST I dD 



2p 



1 + (av) 



NsT(M,t). (16) 



Note that in this case, because of the extra term, the halo formation rate can be negative for some 
values of halo mass. Since negative values of rate of halo formation are unphysical, this indicates 
that the generalization of Sasaki approximation to the Sheth-Tormen mass function is incorrect 



The s ame problem is encountered if we use other models of the halo mass function (iSamui et al 



120091) . 



However, since the basic framework outlined in the beginning of this section is clearly cor- 
rect, there should not be any problems in generalizing it to other mass functions. It is therefore 
likely that the simplifying assumptions of the Sasaki method that led to the estimate of the halo 
destruction rate efficiency of equation (fTsT i are responsible for negative halo formation rate. 



2.3 Excursion set approach to halo formation rates: Press-Schechter mass function 

To check this assertion we perform an explicit calculation of the rate of halo formation using 
the excursion set formalism. Recall that from equations (|2) and (O, we can write for the halo 
destruction rate efficiency as 

oo 

= (.7, 

Ml 

where Q{M[,M2', t) represents the probability that an object of mass Mi grows into an object of 
mass M2 per unit time through merger or accretion at time t. This quantity is also known as the 
transition rate. 



In the excursion set formalism, the conditional probability for a halo of mass Mi present at 
time ti to me rge with another halo to fo rm a larger halo of mass between M2 and M2 + dM2 at 
time t2 > t\ (ILacev & Colelll993l 1 1 9941) can be written for the extended Press-Schechter mass 
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function as 



f{M2,62\Mu5x)dM2 



X exp 



2 ^2(^1 

71 5 

2\2 



■Si) 2 



(520-f - 610-^) 

lalcrlicrl-crl) 



do-2 



dM2 



(18) 



Here, cri and 0-2 are values of the standard deviation of the density perturbations when smoothed 
over scales that contain masses Mi and M2 respectively, and 61 and 62 are the values of the 
threshold density contrast for spherical collapse at time fi and f2 respectively. Taking the limit t2 
tends to fi, i. e. 62 tends to 6], we can determine the mean transition rate at time t - ti: 



Q(Mi,M2;t)dM2 



X exp 



n 



■2 1 ^2 ^2 



lcr\(T\ 



da-2 





d5 




dt 



dM2 



dM2. 



(19) 



This represents the probability that a halo of mass Mi will accrete or merge to form another halo 
of mass M2 at time f. We can use this with equation (fTTb to explicitly compute the destruction 
rate, and hence the halo formation rate. 



However, in the excursion set method, an arbitrarily small change in the halo mass is treated 
as creation of a new halo. As a result, the integral in equation ( flTl i diverges unless we specify 
a "tolerance" parameter. We assume that a halo is assumed to have survived unless its mass in- 
creases such that Ml — > M2 > Mi(l + e) due to either accretion or merging, where e is a small 
number This assumption allows us to introduce a lower cutoff in the integral in equation ( fTTl i 
and the lower hmit changes to Mi(l + e), leading to a convergent integral. This is also physic- 
ally pertinent for our application as infinitesimal changes do not lead to variations in dynamical 
structure of haloes, and hence we do not expect any changes in galaxies hosted in haloes that do 
not undergo a major merger This is similar in spirit to the assumption made el sewhere in the 
literature that a halo is assumed to survive until its mass increases by a factor two (Lacev & Cole 
fl994: Kita yama & Su to 1996). Note that N-body simulations have a natural cutoff due to the 
discrete nature of N-body particles. With the introduction of this new parameter, the modified 
formula for the halo destruction rate efficiency is given by 

0(Mi,f)= J Q{Mi,M2\t)dM2 (20) 

This can then be used to calculate the rate of halo formation using equation Q. 

Fig. [U shows the destruction rate efficiency (piM, t) computed in this manner for the Press- 
Schechter mass function for an Einstein-de Sitter cosmology with power law spectrum of density 
perturbations with index -1.5. Curves have been plotted for e = 0.1 and e = 0.5. We have also 
shown the Sasaki approximation in the same panel. The excursion set result has three features: 
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0.1 



0.01 



z = 10.0 (n = -1.5) 



PS (e = 0.1) 
PS (c = 0.5) 



\ 

\ 

\ 

\ 



10 * 10 ^ 0.01 0.1 1 10 100 1000 10* 10" 

Figure 1. Destruction rate 0(M, f) at z = 10 for the Press-Schechter mass function for a power law model 
with index -1.5. Curves have been plotted for e = 0.1 and e = 0.5. 

1 . At small M, the excursion set value approaches the destruction rate computed using the 
Sasaki approximation. 

2. The destruction rate has a peak, more pronounced for smaller e, near the scale of non- 
linearity. 

3. At larger scales the destruction rate falls rapidly; this is the region where deviations from 
the Sasaki result are the largest. Thus the halo destruction rate efficiency vanishes at large 
masses. 



A similar trend is seen for other power spectra. We postpone a detailed discussion of these issues 
to the end of this section. 



2.4 Excursion set approach to halo formation rates: Sheth-Tormen mass function 



As discussed in Subsection l2.2l the Sheth-Tormen mass function is known to be a much better fit 
to N-body simulations than the Press-Schechter mass function. Several other forms of halo mass 



funct i on have also been fitted to results o f high resolution N-body simulations (IJenkins et al 



20011 : iReed et alj|2003t I Warren et al.ll2006h . In this paper we focus on the Sheth-Tormen mass 



function. Recall that the Sasaki prescription gave unphysical results when applied to this form 
of the mass function. Therefore, we now derive the halo destruction rate efliciency, and the halo 
formation rates for the Sheth-Tormen mass function. This requires obtaining analogs of equations 
(HI and dull. 
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Sheth. Mo & TormeDl(l200l l) showed that once the barrier shape is known, all the predictions 
of the excursion set approach, like the conditional mass function, associated with that barrier can 
be computed easiljH Further, they showed that the barrier shape for ellipsoidal collapse is 



B(cr, t) = 6,Jcr, t) = ^/^6c(t) [1 + /^(av)-^] ., 



(21) 



where a - 0.75, 13 - 0.485, y - 0.615, and, SM) is th e threshold value of overdensity required 
for spherical collapse (also see lSheth & Tonnenll2002l) . They also found that, for various barrier 
shapes B{S ), the first-crossing distribution of the excursion set theory is well approximated by 



f(S)dS 



\T(S)\ 



exp 



B(Sf 



2S 



dS, 



where T(S ) denotes the sum of the first few terms in the Taylor series expansion of B(S ) 



(-S)" d"B(S) 
~. dS" ■ 



(22) 



(23) 



(Here, for conformity with the literature, we use the symbol S = cr^.) This expression gives the 
exact answer in the case of constant and linear barriers. For the ellipsoidal barrier, we can get 
convergence of the numerical result if we retain terms in the Taylor expansion up to « = 5. 



For Press-Schechter mass function, the conditional mass function f(S i,6i\S2, ^2) can be ob- 
tained from the first crossing f{S ) by just changing the variables 6 ^ 61 - 62 and S ^ 5 1 - 5 2. 
This can be done because, despite the shift in the origin, the second barrier is still one of constant 
height. This is no longer true for Ellipsoidal collapse and hence we cannot simply rescale the 
function of equation to get the conditional mass function. Instead, this can be done by making 
the replacements B{S ) — > B(S 1 ) - B(S 2) and 5 —> 5 1 - 5 2 in equation (|22]) . 

\T(Si\S2)\ 



f(Si\S2)dSi 



V27r(5 1-52)3/2 

iBiSi)-BiS2)f 



X exp 



2(5 1-52) 



dS I, 



where we now have 



T(Si\S2)^J] 



{-(Si-S2)r d" (B(5i)-B(52)) 



dS" 



(24) 



(25) 



Using Bayes' theorem, we now have 

\T(Si\S2)\\T(S2) 



f{S2\Sl)dS2 = 



\TiSi)\ 



V2^ 



Si 



X exp 



[B(Si)-B{S2)]- 
2(5i-52) 



52 (5 1-52) 
2, 



BHS2) B'jSi) 
252 ^ 25 1 



dSo 



(26) 



^These can be calculated for non-Gaussian initial conditions, see, e.g., lDe Simone. Maggiore & Riottd j201lh 
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Figure 2. Same as Fig.[T]but for the Sheth-Tormen (as well as Press-Schechter) mass function. 

A change of variables from 5 to M now gives us an analog of equation (fTST l for the Sheth-Tormen 
mass function. In other words, we get the conditional probability fsY{M2\M\)d\vi M2 that a halo 
of mass Ml present at time t\ will merge to form a halo of mass between M2 and M2+dM2 at time 
t2 > ti. Further, taking the limit as t2 tends to ti{= t), we obtain Q{Mi,M2', t). As before, we can 
then use it to calculate the halo destruction rate efficiency 0(M, t) and the rate of halo formation 
Af|^]^(Mi,z) using equations (|5]l and (l20l l. We perform this part of the calculation numerically. 
It is also possible to use this formalism to calculate formation rates for the square-root barrier 
( Moreno. Giocoli & Shethl[2009l I2OO8 : Giocoli et al. 20071) . which is a good approximation for 
the ellipsoidal collapse model. We do not attempt this calculation here as it is beyond the scope 
of this paper 

Fig. |2]is the analog of Fig. [T]for the Sheth-Tormen mass function. It shows the destruction 
rate per halo 0(M, t) computed using the excursion set method for an Einstein-de Sitter cosmology 
with power law spectrum of density perturbations with index -1 .5 at z = 10.0. Curves have been 
plotted for e = 0.1 and e - 0.5. We have also shown the Sasaki approximation for ST mass 
function in the same panel for comparison. This result for the Sheth-Tormen mass function 
has the same features as the result for the Press-Schechter mass function. We also see that the 
destruction rate efficiency is far from constant at small M/M„i. Thus the central assumption of 
the Sasaki prescription is invalid in the case of Sheth-Tormen mass function as well. 



3. N-body simulations 



From the excursion set calculation described in the previous section, we thus find that the halo 
destruction rate efficiency is not independent of mass as is assumed in the Sasaki prescription. 
Clearly, this is the reason why Sasaki prescription yields unphysical values for the rate of halo 
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Table 1. For power law: here n is the index of the power spectrum, A'box is the size of the simulation box , 
A'piirt represents the number of particles, rj^j is the scale of non-linearity at the earliest epoch, r^j is the actual 
scale of non-linearity for the last epoch, r™" represents the maximum scale of non-linearity and z, is the 
starting redshift of the simulations for every model. 



n 


A^box 


A^part 




nl 


'nl 


Zi 


-1.5 


400^ 


400^^ 


2.5 


12.0 


10.0 


103.38 


-0.5 


256^ 


256^^ 


2.5 


12.0 


18.2 


291.53 



Table 2. For LCDM: columns 1 and 2 list the size of the box (in Mpc//i) and the number of particles used 
in the simulations. Columns 3 and 4 give the mass (in Mq/H) and force resolution (in kpc//i; not to be 
confused with the e used in the text) of the simulations, while columns 5 and 6 tell us the redshift at which 
the simulations were terminated and the redshift for which the analyses were done. 



^box 


A^part 


W^pait 


e 


Zf 


^out 


23.04 


512^ 


6.7 X lO*- 


1.35 


5.0 


5.04 


51.20 


512^ 


7x 10^ 


3.00 


3.0 


3.34 


76.80 


512^ 


2.3 X 10*^ 


4.50 


1.0 


1.33 



foiTnation. In this section and the next, we now compare the results of our excursion set calcula- 
tion with results of N-body simulations. 



W e used the TreePM code (lKhandai& Baglall2009h for these simulations. The TreePM 



(lBaglall2002 : iBagla & Ravll2003r) is a hybrid N-body method w hich im proves the accuracy and 

performanc e of the Barnes-Hut (BH) Tree method (Barnes & HutI 198 6|) by combining it with the 

PM inethod ("Miller 1983; Klypin & Shandarin 1983; Bouchet, Adam & Pellat 1985; Bouc het & Kandrupl 
1 19851; iHocknev & Eastwood,1988^ .Bagla & Padmanabhan,1997: ,Merz. Pen & Trac.,2005 i The 
TreePM method explicitly breaks the potential into a short-range and a long-range component 
at a scale r,: the PM method is used to calculate the long-range force and the short-range force 
is computed using the BH Tree method. Use of the BH Tree for short-range force calculation 
enhances the force resolution as compared to the PM method. 

The mean interparticle separation between particles in the simulations used here is Z,nean = 1-0 
in units of the grid-size used for the PM part of the force calculation. In our notation this is also 
cube root of the ratio of simulation volume A^^^^^ to the total number of particles A^pan. 

Power law models do not have any intrinsic scale apart from the scale of non-linearity in- 
troduced by gravity. We can therefore identify an epoch in terms of the scale of non-linearity 
r„\. This is defined as the scale for which the linearly extrapolated value of the mass vari- 
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ance at a given epoch criia, r^x) is unity. All power law simulations are normalized such that 
cr^(fl - 1.0, Tni = 8.0) = 1.0. The softening length in grid units is 0.03 in all runs. 

The ACDM simulations were run with the set of cosmological parameters favored by Wilkin- 
son Microwave Anisotropy Probe 5-yr data (WMAP; Komatsu et al. 2009) as the best fit for the 
ACDM class of models: Q„, = 0.2565, Qa = 0.7435, = 0.963,0-8 = 0.796, /i = 0.719 and 
Qfo/j^ = 0.02273. The simulations were done with 5 12-' particles in a comoving cube of three 
different values of the physical volume as given in Table |2] 

Simulations introduce an inner and an outer scale in the problem and in most cases we work 
with simulation results where Lbox s> r,,] > Z-erid, where Lgdd, the size of a grid cell is the inner 
scale in the problem. Lbox is the size of the simulation and represents the outer scale. In Table 
([TJ we list the power law models simulated for the present study. We list the index of the power 
spectrum n (column 1), size of the simulation box A^box (column 2), number of particles A^pan 
(column 3), the scale of non-Unearity at the earliest epoch used in this study (column 4), and, 
the maximum scale of non-linearity, r™" (column 6) given our tolerance level of 3% error in 
the mass variance at this scale. For some models with very negative indices we have run the 
simulations beyond this epoch. This can be seen in column 5 where we list the actual scale of 
non-linearity for the last epoch. The counts of haloes in low mass bins are relatively unaffected by 
finite box considerations. We therefore limit errors in the mass function by running the simulation 
up to r™^" . Column 7 lists the starting redshift of the simulations for every model. Similarly, in 
Table @, we mention the details of the LCDM simulations used in this work. We list the size of 
the simulation box Lbox in Mpc//i (column 1), number of particles used in the simulations A^pan 
(column 2), mass of the particles Wpan in Mq/Zi (column 3), force resolution e (not to be confused 
with the e used in the text) of the simulations in kpc//i (column 4), the redshift Zf at which the 
simulations were terminated (column 5) and the redshift Zout for which the analyses were done 
(column 6). 

In order to follow the merger history of dark matter haloes in each of these simulations, we 
store the particle position and velocities at different redshifts. A friend-of-friend group finding 
algorithm is used to locate the virialised haloes in each of these slices. We adopt a linking length 
that is 0.2 times the mean inter-particle separation, corresponding to the density of virialised 
haloes. Only groups containing at least 20 particles are included in our halo catalogs. A merger 
tree is then constructed out of the halo catalogs by tracking the evolution of each particle through 
various slices. This lets us identify a halo as it evolves with time through mergers with other 
haloes. We then describe the formation and destruction of haloes in terms of change in number 
of particles between consecutive snapshots of the simulation. When a halo of mass M at redshift 
z turns into a halo of mass M' at z'{< z), then we say that a halo of mass M was destroyed at 
redshift z and a halo of mass M' has formed at z' if M' > M{\ + e). We identify the resolution 
parameter e with that used in our excursion set calculation and experiment with different values 
as described in the next section. 

We find that a tolerence parameter e, similar to the one defined before, also occurs while 
analysing the results of N-body simulations. We identify these two quantities. As we will see in 
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Figure 3. Comparison of the destruction rate efficiencies computed using our method and Sasaki formalism 
for both ST and PS mass function at r„i = 5 grid lengths (top row) and r„; = 8 grid lengths (second row). 
All curves are plotted for power-law model with index n = -0.5. Curves for e = 0.5 are shown in the left 
panel and 6 = 0.1 in the right panel. Points with error bars represent the corresponding results obtained from 
N-body simulations. 



the next section, the formation rate in our model has a dependence on e, which reproduces the 
dependence of the results of N-body simulations on this quantity. Thus, the presence of e in our 
analytical model is crucial in comparing our results with the N-body results. 



4. Results and discussion 



In this section we present the results of a comparison of our calculations presented in Section |2] 
with N-body simulations. We present comparison of the destruction rate efficiency and the rate 
of halo formation and then discuss our results at the end of this section. We also consider two 
related quantities, the halo survival probability and the distribution of halo formation times, that 
were defined in Section|2] 
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Figure 4. Same as Fig. [3] but now for n = -1.5. The two epochs correspond to r„i = 4 and r„i = 8 grid 
lengths respectively. 

4.1 Halo destruction rate efficiency 



Figs. [3] and |4] show the halo destruction rate efficiency (p(M, f) for Sheth-Tormen and Press- 
Schechter mass functions in an Einstein-de Sitter universe with a power law power spectrum 
of density fluctuations with indices n - -0.5 and n = -1.5 respectively. The top row of both 
figures shows the halo destruction rate efficiency at z = 0.8 and the second row shows the same 
at z = 0.0. In each case, we compute the halo destruction rate efficiency using the Sasaki method 
as well as our exc ursion set method. We then derive (piM, f) from our N-body simulations for 
a comparison: see iBagla. Khandai & Kulkarnj ( EoO^ for details of the simulations and best fit 
parameters for the ST mass function. These results are superimposed on the plots. For the excur- 
sion set calculation and for the comparison with simulations, we use e = 0.5 (left column) and 
e = 0.1 (right column). For the two power spectra, the two redshifts that we consider correspond 
to Tni - 5 and r^x - 8 grid lengths, and r,ii = 4 and r,ii = 8 grid lengths respectively. 



As we saw in Figs. [T] and |2l we find that Sasaki's assumption is not valid for ST or PS 
mass functions, that is 0(M, t) depends on the halo mass. We also see that the value of (p{M, t) 
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Figure 5. Destruction rates for ACDM model for both PS and ST mass functions using different thresholds 
(6 = 0.5 for first; e = 0.1 for second) and different redshifts (z = 10.2 for left panel, z = 2.0 for right panel). 
Again, points with error bars represent the corresponding results obtained from N-body simulations. 



derived from simulations matches well with that calculated by our method. On the other hand, 
the predictions of Sasaki's approximation do not match the simulations. This difference is more 
pronounced for the smaller value of e. Note that the points from N-body simulations have large 
error-bars at higher mass as the number of haloes decreases at these scales. The most notable 
feature of the destruction rate efficiency in the excursion set picture is that it cuts off very sharply 
for large masses. Another aspect is that for small e, there is a pronounced peak in (p and it drops 
off towards smaller masses. 

We have also calculated the destruction rate efficiency for the ACDM cosmological model 
for both Press-Schechter and Sheth-Tormen mass functions and compared it with derived values 
from simulations. The results are shown in Fig. |5]for two redshifts (2.0 and 10.2) and two values 
of e (0.5 and 0.1). We can see that results calculated by our technique fit numerical results better. 
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Figure 6. Comparison of the foitnation rates computed using our metliod and Sasald formalism for both ST 
and PS mass functions for r„i = 5 (top row) and r„i = 8 (bottom row). All curves are plotted for power-law 
model with index n = -0.5. Curves for 6 = 0.5 are shown in the left panel and 6 = 0.1 in the right panel. 
Points with error bars represent the coiTesponding results obtained from N-body simulations. 



4.2 Halo formation rate 



Having calculated the destruction rate efficiency, we can now calculate the halo formation rate 
using the formalism described in Section |2] and compare it with the derived halo formation rates 
from our simulations. The results are shown in Figs. |6]and|7]for an Einstein-de Sitter Universe 
with a power law power spectrum of density fluctuations with indices n - -0.5 and n = - 1 .5 
respectively. The third row of both figures shows the formation rate at redshift z = 0.8 and 
the fourth row shows the same at redshift z = 0.0. Note the quantity plotted here is the ratio 
NformiM, t)/N(M, t). We have shown the results from the Sasaki prescription and the excursion 
set calculations and have superimposed formation rates derived from N-body simulations. As 
before, for the excursion set calculation and for the comparison with simulations, we use e = 0.5 
(left column) and e = 0.1 (right column). For the two power spectra, the two redshifts that we 
consider correspond to r^x - 5 and rni - 8 grid lengths, and r^x - 4 and r,,] - 8 grid lengths 
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Figure 7. Same as Fig. |6] but now for n = -1.5. The two epochs correspond to r„i = 4 and r„i = 8 grid 
lengths respectively. 



respectively. The corresponding results for the ACDM cosmological model are shown in Fig. [8] 
for two redshifts (2.0 and 10.2) and two values of e (0.5 and 0.1). 

Again, we see that the excursion set results fit simulation data much better as compared to 
the results from Sasaki prescription. The Sasaki method underestimates the formation rates by a 
large factor for low mass haloes. Results from the two methods tend to converge in the large mass 
limit, although a systematic difference remains between the Sheth-Tormen and Press-Schechter 
estimates, with the former always being larger that the later The difference in the Sasaki estimate 
and the excursion set estimate for the destruction rate efficiency and the formation rate is as 
high as an order of magnitude at some scales so the close proximity of simulation points to the 
excursion set calculations is a clear vindication of our approach. It is worth noting that there 
is a clear deviation of simulation points from the theoretical curves at small mass scales and 
this deviation is more pronounced at small mass scales for e = 0.5. It may be that some of the 
deviations arise due to a series representation of the barrier shape, and the number of terms taken 
into account may not suffice for the estimate. We have found that truncation of the series can 
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Figure 8. Formation rates for ACDM model for both PS and ST mass functions using different thresholds 
(6 = 0.5 for first; 6 = 0. 1 for second row) and different redshifts (z = 10.2 for left panel, z = 2.0 for 
right panel). As usual, points with error bars represent the corresponding results obtained from N-body 
simulations. 



affect results at small masses, though in most cases results converge with the five terms that we 
have taken into account for the range of masses considered here. 



4.3 Halo survival probability 



An important auxiliary quantity in the ongoing discussion is the halo survival probability, defined 
in Section |2] From our calculation of the halo destruction rate efficiency, we calculated the 
survival probability of dark matter haloes using both the excursion set formalism and the Sasaki 
prescription and compared results. These results are shown in Fig. |9] which shows the survival 
probabilities in the ACDM cosmological model for the Press-Schechter (left panel) and Sheth- 
Tormen (right panel) mass functions using the two approaches at three different redshifts (z = 1, 
3 and 10). In this case, we have used e = 0.5 for the excursion set calculation. 
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Figure 9. Comparison of the survival probabilities computed using our method and Sasaki formalism for 
both PS (left panel) and ST (right panel) mass functions with different redshifts (z = 1,3 and 10) for e = 0.5. 
Curves have been plotted for the ACDM model. These curves show the probability that the halo survives 
from that redshift up to the present epoch. 



In Sasaki approximation, the destruction rate is independent of mass and hence the survival 
probabiHty is also independent of mass. Our calculations show that this approximation is not true, 
and hence the survival probability of haloes must also depend on mass. We note that the survival 
probability is high for large mass haloes: if a very large mass halo forms at a high redshift then 
it is likely to survive without a significant addition to its mass. Smaller haloes are highly likely 
to merge or accrete enough mass and hence do not survive for long periods. Survival probability 
drops very rapidly as we go to smaller masses. While this is expected on physical grounds, it is 
an aspect not captured by the Sasaki approximation where equal survival probability is assigned 
to haloes of all masses. The mass depen dence of survival probability is qualitatively similar to 
that obtained bv lKitavama & Sutol ( ll99d) . There is no significant qualitative difference between 
the curves for the Press-Schechter and the Sheth-Tormen mass functions. 



4.4 Formation time distribution 



Finally, another interesting quantity is the distribution F(M; tf, t) of formation epochs tf of haloes 
with mass M at f, defined in Section |2] This distribution can be obtained once the survival 
probability and formation rate of haloes is known. We calculated the formation time distribution 
using the excursion set formalism and the Sasaki prescription. The results are shown in Fig.fTOl 
We plot F(M;zf,z - Q)/N(M,z - 0) versus the formation redshift Zf for three different masses 
(10'^ 10'" and IO^Mq) in the standard ACDM model for both Press-Schechter (left column) and 
Sheth-Tormen (right column) mass functions with e - 0.5 (first row) and e = 0.1 (second row). 
A common feature is that F as a function of Zf increases up to a certain redshift and then starts to 
decline. The epochs at which F drops by an order of magnitude from its peak can be interpreted 
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Figure 10. Plots for formation epoch distribution of lialoes. Left column is for the PS and the right column is 
for the ST mass function. Curves have been plotted for the LCDM model. The formation epoch distribution 
as computed using the Sasaki formalism and the excursion set approach described in this work is shown in 
the top panel. 



as typical range of redshifts for the formation of bound systems of respective masses which exist 
at z = 0. 

The differences between the formation redshift distribution for e - 0.5 and e = 0.1 are along 
expected lines: the formation redshifts are smaller for the lower value of e as a smaller change 
in mass is required for us to declare that a new halo has formed and hence typical haloes do not 
survive for a very long time. We see that the excursion set calculation suggests that haloes formed 
more recently as compared to the Sasaki approximation based estimate. This can be understood 
in terms of the equal survival probability assigned by the Sasaki approximation to haloes of all 
masses. For a clearer comparison, the ratio of the estimate based on Sasaki approximation and 
the excursion set calculation is shown in Fig.[TT] We note that for very low mass haloes these two 
estimates differ by more than an order of magnitude. The main qualitative difference between 
the plots for the Press-Schechter and the Sheth-Tormen mass functions is caused by the negative 
formation rates in the Sasaki approximation. 



Halo formation rates 



23 




4.5 Discussion 



The results described above show conclusively that the excursion set approach predicts halo form- 
ation and destruction rates that match with simulations much better than the Sasaki approxima- 
tion. 

Another noteworthy aspect is that the destruction and formation rates depend on the value of e 
in simulations as well as the excursion set calculation thereby allowing us to differentiate between 
major and minor mergers. In comparison, there is no natural way to bring in this dependence in 
the Sasaki approximation. While the match between simulations and the excursion set approach 
for the two values of e is satisfying, it raises the question of the appropriate value of this parameter 
In our view the appropriate value of the parameter should depend on the application in hand. 
In semi-analytic galaxy formation models, we should use a value of e that corresponds to the 
smallest ratio of masses of the infalling galaxy and the hos t galaxy where we expec t a significant 



dynamical influence on star formation rate. For instance, iKauffmann et al.l (119991) use e = 0.3 
in their semi-analytic galaxy formation model while considering formation of bulges in merger 
remnants. In case of galaxy clusters we may base this on the smallest ratio of masses where the 
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Figure 12. Upper panels show foiination rates for ST mass function. Lower panels show the same where 
we used (p computed from excursion set approach in the PS mass function and used that to compute the 
formation rate in the ST mass function. 
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Figure 13. Ratio of formation rates estimated in the two approaches shown in Fig.ll2l 
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Figure 14. Same as Fig.[T2]but for 6 = 0.1 




Figure 15. Same as Fig.[T3]but for e = 0.1 
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intra-cluster medium is likely to be distur bed in a manner acce ssible to obse rvations in X-ray 



emission or the Sunyaev-Zel'dovich effect (ISunvaev & Zeldovich> 1972; .Navarro. Frenk & White 
ll995l:iKavi[2004h . 



While the close match between simulations and the excursion set calculation is useful, it also 
implies that we should not use the simpler Sasaki approximation. The excursion set calculation 
of the halo destruction rate is fairly simple for the Press-Schechter mass function, but the cor- 
responding calculation for the Sheth-Tormen mass function is much more complicated. Plots of 
the destruction rate efficiency (p{M) for all the models suggest that its variation with mass and e 
is very similar for the PS and ST mass function. This suggests an approximation where we use 
4>{M, z; e) computed using the Press-Schechter mass function and use that to compute the halo 
formation rate in the Sheth-Tormen mass function. Figs. [12] and [14] show the halo formation rate 
for the ACDM model at different redshifts and compare the excursion set calculation, the Sasaki 
approximation and the intermediate approximation suggested above. We have also shown the 
ratios of formation rates estimated in these two approaches mentioned above in Figs. [T3]and[T5] 
We find that the intermediate approximation is not plagued by negative halo formation rates and 
that it is an excellent approximation at all mass scales at higher redshifts. At lower redshifts, the 
approximation is still good at high masses but not so at smaller masses. 

While comparing our analytical results with those of N-body simulations, we find a system- 
atic deviation between the two at the high mass end. This is possibly related to the problm of 
'halo fragmentation' while deriving halo merger trees from the simulations. In about 5% of all 
haloes, particles in a given progenitor halo can become part of two independent haloes at a future 
epoch. This is usually attributed to the fact that the FOF algorithm groups particles based on 
the inter-particle distance. This can result in the identification of two haloes separated by a thin 
'bridge' of particles to be treated as a single halo. Such hal o fragmentation has bee n treated using 



different techniques in various halo formation rate studies. iFakhouri & Ma (12008 ') compare these 



techniques and find that the effect of halo fragmentation is maximum of high mass haloes. 



5. Conclusions 



Key points presented in this paper can be summarized as follows: 



• We revisit the Sasaki approximation for computing the halo formation rate and compute 
the destruction rate explicitly using the excursion set approach. 

• We introduce a parameter e, the smallest fractional change in mass of a halo before we 
consider it as destruction of the old halo and formation of a new halo. 

• We show that the halo destruction rate is not independent of mass even for power law 
models and hence the basis for the Sasaki ansatz does not hold. Two prominent features of 
the halo destruction rate are the rapid fall at large masses, and a pronounced peak close to 
the scale of non-linearity. The peak is more pronounced for smaller values of e. 
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• Using the excursion set approach for the Sheth-Tormen mass function leads to positive 
halo formation rates, unlike the generalization of the Sasaki ansatz where formation rates 
at some mass scales are negative. 

• We compare the destruction rate and the halo formation rates computed using the excur- 
sion set approach with N-body simulations. We find that our approach matches well with 
simulations for all models, at all redshifts and also for different values of e. 

• In some cases there are deviations between the simulations and the theoretical estimate. 
However, these deviations are much smaller for the excursion set based method as com- 
pared to the Sasaki estimate. 

• It may be that some of the deviations arise due to a series representation of the barrier 
shape, and the number of terms taken into account may not suffice for the estimate. We 
have found that truncation of the series can affect results at small masses, though in most 
cases results converge with the five terms that we have taken into account for the range of 
masses considered here. 

• We show that we can use the halo destruction rate computed for the Press-Schechter mass 
function to make an approximate estimate of the halo formation rate in Sheth-Tormen mass 
function using equation (|5]l. This approximate estimate is fairly accurate at all mass scales 
in the ACDM model at high redshifts. 

• The halo survival probability is a strong function of mass of haloes, unlike the mass inde- 
pendent survival probability obtained in the Sasaki approximation. 

• The halo formation redshift distribution for haloes of different masses is also very different 
from that obtained using the Sasaki approximation. This is especially true for the Sheth- 
Tormen mass function where the Sasaki approximation gives negative halo formation rates 
in some range of mass scales and redshifts. 



The formalism used here for calculation of halo formation rate and other related quantities 
can be generalized to any description of the mass function if the relevant probabilities can be 
calculated. Within the framework of the universal approach to mass functions, it can also be 
used to study form ation rates of haloes in different cosmological models dLinder & Jenkinsll2003l 



Maccio et al This allows for an easy co mparison of theory with obse rvations for quantities 



like the major merger rate for galaxy clusters (ICohn. Bagla & Whitell2001 ) 



In case of semi-analytic models of galaxy formation, our approach allows for a nuanced 
treatment where every merger need not be treated as a major merger and we may only consider 
instances where mass ratios are larger than a critical value for any affect on star formation in the 
central galaxy. 
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